
function outcomes=inference(turnouts, demographics, propensities)

    % 2011.9.13: The code is based on ATE_TT.m 
    % I rewrite in along with eco_bootstrap.m
    
    
    %% Input
    % turnouts: stateicp, year(binary), republican votes, total voters.
    % demographics: stateicp, year(binary), log income, education, age, white, unemployed.
    % propensities: propensity score of year==1.
    % Note that each rows in 'demographis' and 'propensities' correspond to one individual sample.
    
    %% Output
    % outcomes=[mu0L, mu0U, mu1L, mu1U, mu01L, mu01U,mu10L, mu10U, YforXD(:,1:2)];
    
    
    % Process Data
    yxd=[(turnouts(:,3)./turnouts(:,4)) , turnouts(:,[1,2])];
    % yxd= rep_vote share(percentage), state(or national), year
    xzd=[demographics(:,1:2), 1./propensities(:,1), 1./(1-propensities(:,1)) , propensities(:,1)./(1-propensities(:,1)) ] ;
    xzd_national=[demographics(:,1:2), 1./propensities(:,2), 1./(1-propensities(:,2)) , propensities(:,2)./(1-propensities(:,2)) ];
    % xzd= state, year(dummy), W, V, VoverW w.r.t propensity scores conditioned on stateicp
    % xzd_national= state, year(dummy), W, V, VoverW w.r.t propensity scores unconditioned on stateicp
    

    % Caution! "yxd" has an additional state with icp number 99. This is for nation-wide voting data.
    
    % For each state and nation-wide(X) and year(D), store values of W,V,VoverW. 
    % We will use this later for computing FHatW, FHatV, and FHatVoverW.

    states=unique(yxd(:,2)); % Caution! states also include nation-wide votes.
    NumofState=length(states); %Caution! NumbofState = the number of states in US + 1.

    global WforXD VforXD VWforXD WVforXD
    YforXD=zeros(NumofState,2);
    WforXD=cell(NumofState,2);
    VforXD=cell(NumofState,2);
    VWforXD=cell(NumofState,2);

    for x=1:NumofState
        if x < NumofState
            for j=1:2 % j is dummy(D) +1.
                
                ind1=find((yxd(:,2)==states(x))&(yxd(:,3)==j-1));
                YforXD(x,j)=yxd(ind1,1);


                ind2=find(xzd(:,1)==states(x)& xzd(:,2)==j-1);

                WforXD{x,j}=sort(xzd(ind2,3));
                VforXD{x,j}=sort(xzd(ind2,4));
                VWforXD{x,j}=sort(xzd(ind2,5));
                WVforXD{x,j}=sort(1./xzd(ind2,5));

            end
        else
            for j=1:2 % j is dummy(D) +1.


                ind1=find((yxd(:,2)==states(x))&(yxd(:,3)==j-1));
                YforXD(x,j)=yxd(ind1,1);


                ind2=find(xzd_national(:,2)==j-1);

                WforXD{x,j}=sort(xzd_national(ind2,3));
                VforXD{x,j}=sort(xzd_national(ind2,4));
                VWforXD{x,j}=sort(xzd_national(ind2,5));
            end
        end

    end



    % Compute Bounds. We first calculate summations over each state(X).
    % Later, we divide them by the number of total observations (m)
    % to obtain unconditional ATE, etc.

    mu1L=zeros(NumofState,1); 
    mu1U=zeros(NumofState,1); 
    mu0L=zeros(NumofState,1); 
    mu0U=zeros(NumofState,1); 
    mu01L=zeros(NumofState,1); 
    mu01U=zeros(NumofState,1);
    mu10L=zeros(NumofState,1); 
    mu10U=zeros(NumofState,1);

    SampleSize=cellfun('length',WforXD);
    

    for x=1:NumofState
        
        % if D=1, we only update mu1L, mu1U

        j=2;
        % We define the functions in the integral.
        mu1L(x) = quad(@(u) ((1-u) > 1-YforXD(x,j)).*InvFHatW(u,x,j) , 0.05, 0.95)*SampleSize(x,2)/sum(SampleSize(x,:));
        mu1U(x) = quad(@(u) (u > 1-YforXD(x,j)).*InvFHatW(u,x,j), 0.05, 0.95)*SampleSize(x,2)/sum(SampleSize(x,:));
        mu10L(x) = quad(@(u) ((1-u) > 1-YforXD(x,j)).*InvFHatVW(u,x,j) , 0.05, 0.95).*(SampleSize(x,1)./SampleSize(x,2));
        mu10U(x) = quad(@(u) (u > 1-YforXD(x,j)).*InvFHatVW(u,x,j) , 0.05, 0.95).*(SampleSize(x,1)./SampleSize(x,2));


        % if D=0, we only update mu0L, mu0U, mu01L, and mu01U.

        j=1;

        mu0L(x) = quad(@(u) ((1-u) > 1-YforXD(x,j)).*InvFHatV(u,x,j) ,0.05,0.95)*SampleSize(x,1)/sum(SampleSize(x,:));
        mu0U(x) = quad(@(u) (u > 1-YforXD(x,j)).*InvFHatV(u,x,j) ,0.05,0.95)*SampleSize(x,1)/sum(SampleSize(x,:));
        mu01L(x) = quad(@(u) ((1-u) > 1-YforXD(x,j)).*InvFHatVW(u,x,j) , 0.05, 0.95).*(SampleSize(x,1)./SampleSize(x,2));
        mu01U(x) = quad(@(u) (u > 1-YforXD(x,j)).*InvFHatVW(u,x,j) , 0.05, 0.95).*(SampleSize(x,1)./SampleSize(x,2));


    end

    outcomes=[mu0L, mu0U, mu1L, mu1U, mu01L, mu01U, mu10L, mu10U, YforXD(:,1:2)];

    
    
end
